import pandas as pd
from scipy import stats
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np


def get_h5_data():
    """
    读取h5的数据
    :return:
    """
    h5 = pd.HDFStore('./data/vstoxx_data_31032014.h5')
    futures_data = h5['futures_data']
    options_data = h5['options_data']
    h5.close()

    return futures_data, options_data


def black_schole_merton(S, r, sigma, X, T):
    """
    bs公式模型
    :param S:
    :param r:
    :param sigma:
    :param X:
    :param T:
    :return:
    """
    S = float(S)
    log_soverx = np.log(S / X)
    half_sigmal_square = 0.5 * np.power(sigma, 2)
    sigma_sqrt_T = sigma * np.sqrt(T)

    d1 = log_soverx + ((r + half_sigmal_square) * T) / sigma_sqrt_T
    d2 = log_soverx + ((r - half_sigmal_square) * T) / sigma_sqrt_T

    return S * stats.norm.cdf(d1, 0.0, 1.0) - X * np.exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0)


def vega(S, r, sigma, X, T):
    """

    :param S:
    :param r:
    :param sigma:
    :param X:
    :param T:
    :return:
    """
    S = float(S)
    log_soverx = np.log(S / X)
    half_sigmal_square = 0.5 * np.power(sigma, 2)
    sigma_sqrt_T = sigma * np.sqrt(T)

    d1 = log_soverx + ((r + half_sigmal_square) * T) / sigma_sqrt_T
    return S * stats.norm.cdf(d1, 0.0, 1.0) * np.sqrt(T)


def impliedvolatility(S, r, sigma_est, X, T, Cstar, it):
    """

    :param S:
    :param r:
    :param sigma_est:
    :param X:
    :param T:
    :param Cstar:
    :param it:
    :return:
    """
    for _ in range(it):
        numer = (black_schole_merton(S, r, sigma_est, X, T) - Cstar)
        denom = vega(S, r, sigma_est, X, T)
        sigma_est -= numer / denom

    return sigma_est


def solve_bs():
    """
    bs期权公式求解
    方法：牛顿迭代法（牛顿-拉夫逊方法）
    """
    # matplotlib中文显示方块
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题

    futures_data, options_data = get_h5_data()
    options_data['IMP_VOL'] = 0.0
    V0 = 17.6639
    r = 0.04
    sigma_est = 2
    tol = 0.5

    for option in options_data.index:
        futureval = futures_data[futures_data['MATURITY'] == options_data.loc[option]['MATURITY']]['PRICE'].values[0]

        if futureval * (1 - tol) < options_data.loc[option]['STRIKE'] < futureval * (1 + tol):
            implied_vol = impliedvolatility(V0, r, sigma_est, options_data.loc[option]['STRIKE'],
                                            options_data.loc[option]['TTM'],
                                            options_data.loc[option]['PRICE'],
                                            it=100)
            options_data['IMP_VOL'].loc[option] = implied_vol

    plot_data = options_data[options_data['IMP_VOL'] > 0]
    maturities = sorted(set(options_data['MATURITY']))

    plt.figure(figsize=(15, 10))

    for maturity in maturities:
        data = plot_data[options_data['MATURITY'] == maturity]

        plot_args = dict(lw=3, markersize=9)
        plt.plot(data['STRIKE'], data['IMP_VOL'], label=pd.to_datetime(maturity),marker='o', **plot_args)

    plt.grid(True)
    plt.xlabel('Strike rate $X$', fontsize=18)
    plt.ylabel('$\sigma$ 的隐含波动率', fontsize=18)
    plt.legend()
    plt.show()


# bs期权公式求解
solve_bs()

